NUMERICAL SIMULATION OF OPTIMAL TRANSPORT PATHS 

QINGLAN XIA 

QQ Abstract. This article provides numerical simulation of an optimal transport 

{**") path from a single source to an atomic measure of equal total mass. We first 

f— ^ construct an initial transport path, and then modify the path as much as 

^vj possible by using both local and global minimization algorithms. 



CO 1. Introduction 

<N 

This article aims at providing numerical simulations of optimal transport paths 
studied earlier in |5],[B],|7] [5] etc. The theory of optimal transport paths was mo- 
tivated by studying a phenomenon in optimal transportation where a transport 



u 

system with a branching structure may be more cost efficient than the one with 
,-h a "linear" structure. Trees, railways, lightning, electric power supply, the circu- 

latory system, the river channel networks, and cardiovascular systems are some 
_h common examples. We used the concept of optimal transport paths between prob- 

ability measures to model such transport systems in |5]. Some related works are 
P] ■ |2] . [3] . [I] and others. Also, in [9], we showed that an optimal transport path is 

t— I exactly a geodesic in the sense of metric geometry on the metric space of probabil- 

J> ity measures with a suitable metric. As a result, people are interested in knowing 

what an optimal transport path look like numerically. This is the motivation of this 
article. Currently, we are using optimal transport paths generated here to model 

X<-\ blood vessel structures found in placentas of human babies and also river channel 

networks. More application of optimal transport paths is expected in modeling 

^— s systems simulating what we seen in the nature. 

qq This article is organized as follows. After briefly recalling some preliminary defi- 

nitions about optimal transport paths, we study algorithms for generating optimal 
transport paths. The idea is as follows. We first consider how to generate an initial 
transport path, then study how to reduce the cost by modifying the initial trans- 
port path as much as possible. We not only use an algorithm of local minimization 






> 

X 

2 but also a global one. Topology of transport paths are not preserved during the 

modification process. These algorithms do not necessarily provide us a perfect opti- 
mal transport path, but it approximate an optimal transport path very well. Using 
eyes of a human being, the author cannot observe a better transport path than the 
path generated here. 

2. Preliminaries 

We first recall some concepts about optimal transport paths between measures 
as studied in [5]. Let X be a convex compact subset in M. d . For any x E X, let 5 X 
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be the Dirac measure centered at x. An atomic measure in X is in the form of 

k 
^ aiS Xi 

with distinct points Xi € X, and Oj > for each i = 1, • • • , fc. Let a and b be two 
fixed atomic measures in the form of 

k i 

(2.1) a = 2_] a>i&xi an d b = V^ bj5 y . 

of equal total mass 

k i 

Definition 2.1. ^4 transport path from a to b is a weighted directed graph G consists 
of a vertex set V (G), a directed edge set E (G) and a weight function 

w :E{G)^> (0,+oo) 

smc/i iftai {#1,2:2,... Xk} U {?/i, j/2, • • • 7 ?/;} C V (G) and /or any vertex v £ V (G) , 

{a i? if v = Xi for some i — 1, • • • , fe 
-fy, ifv = yjforsomej = l,---,l 
0, otherwise 

e~ — v e + —v 

where e~ and e + denotes the starting and ending endpoints of each directed edge 
e€E(G). 



Remark 2.2. The balance equation (2.2) simply means that the total mass flows 
into v equals to the total mass flows out of v. When G is viewed as a polyhedral 
chain or current, (2.2) can be simply expressed as 

<9G = b-a. 

Let Path(&, b) be the space of all transport paths from a to b. 

Definition 2.3. For any a < 1, and any G € Path(&,h), define 

M a (G):= J2 [w(e)] a length (e). 

e£E(G) 

In 5 Proposition 2.1], we showed that for any transport path G G Path (a, b), 
there exists another transport path G G Path (a, b) such that 

M Q (g) < M a (G) , 

vertices ^(G) C V (G) and G contains no cycles. Here, a weighted directed 

graph G — {V (G) , E (G) ,w : E (G) — v (0, 1]} contains a cycle if for some k > 3, 
there exists a list of distinct vertices {i>i,i>2, • • • ,Vk} in V (G) such that for each 
i = 1, • • • , k, either the segment [vi, Vi + {\ or [u,+i, i>i] is a directed edge in E{G), 
with the agreement that v^+i — v\. When a directed graph G contains no cycles, 
it becomes a directed tree. 

An M„ minimizer in Path(&, h) is called an optimal transport path from a, b. 
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3. Simulation of optimal transport paths from a single source 
Let X be a convex subset in W. d . Given two atomic measures in the form of 



A 



(3.1) 



a = mSo and b = > mi5 Vi with to 



i=i 



N 
i=l 



in X of equal total mass, we are interested in seeing what an optimal transport 
path G from the single source a to b look like numerically. 

If N = 1, then G is clearly consisting of only one edge [O, yi) with weight to. If 
N = 2, then we can calculate the optimal transport path as follows. 

3.1. One source to two targets. Suppose there are two atomic measures 

(3.2) /i = mpSo and v = nipSp + toq«5q, 

with too = top + toq for three points P,Q,0 in the space X. To find an optimal 
transport path from /i to v, we need to minimize the function 

f{B) = {m ) a 



OB 



+ (to p )« BP + (m Q ) a 



BQ 



among all points B in the triangle APOQ. Here, we use the notation BP etc. to 

— > 
denote the vector P — B in X, and let BP be the magnitude of this vector. Since 

/ (B) is a continuous function on a compact set, / must achieve its minimum at 
some point B* . Indeed, we can find B* as follows. Suppose B* is located in the 
interior of the triangle APOQ, then it must satisfy the balance equation 



(m ) a 



OB 



OB 



(to p )" -^^- + (m Q ) a -^^- = 



BP 



BQ 



at B = B* . From it, one can easily find the angles 

ZOB*P = 9 1 ,ZOB*Q = 6 2 and IPB*Q = <? 3 



where 
(3.3) 



cos 



k 2 - k x - 1 
2^/kx 



cos 



fci - k 2 - 1 
2^/ko 



cos 



I — ki — k-2 
2 v / fcifc 2 



for 



fcr= — ,fc 2 = — • 
\m J \m J 

Let M (and H) be the projection of the point Q (and P, respectively) along the 

segment OP (and OQ respectively). Then, the centers R (, and £) of the circles 

passing through the triangles AOB*P (and AOB*Q respectively) is given by 



R = 



S 



O + P cot 6>i QM 

2 2 QM 



OP 



O + Q cot 2 P# 
2 2 p~H 



OQ 



QINGLAN XIA 



where 



()M = OPO ® OP - OQ and PH = ° P '°® 0Q - OP. 



OP 



OQ 



By (3.3), 



cot d\ 



fci -1 



and cot 9? 



fci 



1 



'4fci - (fe - *i - 1) V 4fc2 - ( fcl _ fe - l > 

Now, i3* is just the reflection of the point O along the segment RS. That is, 



B* = 2 [(1 - X)R +XS}-0 with A = R ° ' R f 



RS 

whenever B* is located in the interior of the triangle APOQ. Note that a necessary 
condition for B* being located in the interior of the triangle APOQ is the angles 
must satisfy 

ZOQP < $ x ,lOPQ < 2 and IPOQ < 6 3 . 

In case the condition fails, we have three degenerate cases. If the angle /LPOQ > 6* 3 , 
then take B* to be O and we get a "V-shaped" path. If the angle LOQP > &\ and 
IPOQ < 9 3 , then take B* to be Q. If the angle IOPQ > 6 2 and IPOQ < 3 , 
then take B* to be P. 






0.5 1 



0.5 1 



FIGURE 1. Cases for transporting [i to v 



As a result, given /i and v in (3.2 1, we achieved a formula for finding B* . The 



optimal transport path G from [itov has at most three edges: [£?*, P] with weight 
mp, [B*,Q] with weight toq and [0,S*] with weight mo- 
We denote the point B* by V (/x, v) and let 



(3.4) 



g(n,v)=f(0)-f(B*) 



which gives the advantage of taking a "Y-shaped" path over taking a "V-shaped" 
path. 
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3.2. The construction of an initial transport path. When N > 3, it is not 

that easy (or may be impossible) to find an exact solution of an optimal transport 
path. Instead, we would like to find an approximately optimal transport path. The 
idea is to construct an initial transport path G £ Path (a, b) and then modify G 
as much as possible until we can not reduce the cost of G any further. 

3.2.1. The method of transporting small number of points. If N < 2, then we have 
found the optimal transport path G as above. If 2 < N < K for a small given 
number K, then for any pair 1 < i < j < N, let 

cjij = g((m t + mj) So,miS Vi + m 3 (5 9j ) 



where the function g is defined as in (3.4 1. Suppose the maximum of {gr,} is 
achieved at 1 < i* < j* < N. Then, the desired path G is given recursively by 

G=G + mi , [B*, yi ,] + m r [B*,y r ], 



where B* = V ((rrij* + rrij*) So, rni*5 y . t + mj*5 y . t ) is the point in X given by (3.4 1, 
and G is the path from a to b = b — mi*5 Vi „ — mj*5 y .« + (mj* + rrij*) 5b* achieved 
by recursively applying this algorithm. 

3.2.2. The subdivision method. To construct an initial transport path in Path (a, b), 
one may simply take a trivial transport path 



E 



■;??■ 



[0,Vi. 



This is an allowable transport path in Path (a, b). Nevertheless, the degree of the 
vertex O (i.e. the total number of edges in G having O as an endpoint) is N. Then, 
it might become time consuming later for modifying the path at the vertex O when 
N is very large. Instead, we use the following subdivision method to construct an 
initial transport path G s d (a, b), which contains no cycles and has degree at most 
K at every vertex for some given K defined below. 
Let 

"3, if d = 2 



A " 2, if d > 3 

and K = X d where d is the dimension of the ambient space M. d . 
Algorithm (subdivision method): 



Input: two atomic measures a, b in the form of (3.1) and a parameter a < 1; 

Output: a transport path G £ Path (a, b) with degree(v) < K for each v £ V(G). 

If N < K , then we use the method of transporting small number of points 
described above to construct a transport path from a to b. 

If N > K, then let Q be a cube in R d that contains the supports of both a 
and b. We may split the cube Q into totally K — X d smaller cubes {Qi} i=1 of 
size equal to t of the size of Q. For each i = 1, • • ■ , K, let Gi be the path 
G s d (b(Qi) ^c(Qi)i D LQi) from the center c(Qi) of the smaller cube Qi to the re- 
striction of b in Qi achieved by recursively applying this algorithm. Also, let Go be 
the path from a to X)i=i D (Qi) $c, -, by using the method of transporting small 

{Qi) 

number of points. Then, 

K 



G = Y J G l 



i=0 
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provides the desired path G sc i (a, b) from a to b. 



1 




V ** * t ♦♦** ♦* # 


1 


o 



1 2 3 

(a) 






(d) 



Figure 2. (a) A single source (2,-1) and the targeting measure 
represented by 100 random points, (b) An initial transport path 
constructed by the subdivision method. (c)A modified transport 
path achieved by repeated modifying the initial path in (b) using 
the local minimization. (d)An optimal transport path achieved by 
modifying the path in (c) using the global minimization method. 



3.3. Modification of an existing transport path. Now, suppose G is an ex- 
isting transport path from a to b that contains no cycles. We want to modify G 
to reduce the transport cost as much as possible. Before describing algorithms, we 
first introduce some concepts about vertices of an transport path G. 

For any two vertices v , u € V (G), we say that v is an ancestor of u and u is a 
descendant of v, if there exists a list of vertices Vi = v, V2, • • • , v/i-i, Vh = u such 
that each [vi, u»+i] is a directed edge in E (G) for i = 1, • • • , h — 1. Also, if [v, u] is 
a directed edge in E (G), then we say that v is a parent of u and u is a child of v. 
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For each vertex u £ V (G) \ {O}, u has exactly one parent p (u) £ V (G) because 
G contains no cycles and has a single source a = mdo ■ Let to (it) be the associated 
weight on the directed edge \p(u) ,u] in E (G) for each u £ V(G) \ {O}, and also 
set to (O) = to. Note that to (v) > to (u) whenever v is an ancestor of u. Moreover, 
the vertex O is always an ancestor of each u. That is, there exists a list of vertices 
v\,v-2,- ■• ,Vk hi V (G) such that [«i,Vi + i] £ E{G) with v\ = O and Vk = u. Then, 
for each t £ [— m (u) , to (u)], we consider the path 

fc-i 
R (G; t, u) :=G-y^t [vi,v i+1 ] £ Path (a - tS + tS u , b) . 

i=\ 

When t > 0, we say that a mass of t is removed from the path G at vertex u, When 
t < 0, we say that a mass of £ is added to the path G at vertex m. Moreover, the 
potential function of G at a vertex u £ V (G) is defined by 

Pg (p (u) , t) + \p (u) — u\ [to (u) a — (to (u) - t) a ] , u^ O 



(3-5) P G («,*)-, Q) u = 

for t e [— m (u) , to (u)]. Note that Pg (u,t) has the same sign as t. 

3.3.1. local minimization. We first use a local minimization method to modify any 
existing transport path G containing no cycles. 

Input: a transport path G £ Path (a, b) containing no cycles and a < 1; 

Output: a locally optimized path G £ Path (a, b) with M a I G) < M a (G). 

Idea: For each vertex u in G, replace G id (u) by G new (u) whenever M a (G id (u)) > 
M a (G new (u)). 

Here, for each vertex u of G, two transport paths G id (u) and G new (u) are 
defined as follows. Let 

Mc = J^ rn (h) S h and fi P = m (w) <5 p(u ) 

frev(G),p(/i)=« 
be two atomic measures corresponding to the children and the parent of u. Then, 

G id{u) = ^2 m(h)[u,h] + m(u)[p(u) ,u] £ Path([i P ,fi c ) 

heV(G),p(h)=u 

is the union of all weighted edges in G sharing u as their common endpoint. On 
the other hand, one may generate another path G new (u) £ Path (p.p, fie) by using 



the method of transporting small number of points stated in 3.2. 1| 
If 

M a (G i d (u)) > M a {G new (u)) , 
then by replacing G m (u) by G new (u) in G, we get a new path 
G = G- G id (u) + G new (u) £ Path (a, b) 

and M a (g) < M a (G) - M a (G old (u)) + M a (G new («)) < M a (G). So, G is a 

transport path with less cost. Replace G by this modified path G, and continue 
this process for all vertices of G until one can not reduce the cost any further. 

The main drawback of this algorithm is that the result is only local minimization 
rather than global minimization. For instance, edges may intersect with each other. 
Sometimes, using eyes of a human being, one can easily observe a better transport 
path. To overcome these drawbacks, we adopt the following algorithm. 
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3.3.2. global minimization. Now, we introduce the following algorithm of global 
minimization: 



Input: two probability measures a, b in the form of (3.1) and a parameter a < 1; 

Output: an approximately M a optimal transport path G € Path (a, h). 

step 1: construct a transport path G from a to b using the subdivision method; 

step 2: modify the existing path G using the local minimization method; 

step 3: subdivide long edges of G into shorter edges; 

step 4: for each vertex u of G, remove a mass of m (u) at vertex u from the path 
G; change the parent p(u) of u to a better one if possible and then add back a mass 
of m (u) at vertex u. More precisely, 

substep 1 : A list of potential parents of u is defined as 

L (u) = {v *E V (G) : \v — u\ < er, and v is not a descendant of u} , 



where 



P G (u,m(u)) 

a — 



[m (u)} a 
and Pq is defined in (|3.5[). Note that the parent p(u) is automatically in L (u) 



because 

P G (p(u),m(u)) + \p(u) -u\m(u) a 

a = = — — r„ > \p(u) - u\ . 

[m(u)\ 

substep 2: By removing a mass of m (u) at vertex u from the path G, we get 
another path 

G = i? (G; m 0) , u) . 
substep 3: For each »el (u) \ {p(u)}, let 

c(v) = -P G (v,-m(u)), 



where P^, is defined as in (3.5 1 with G replaced by G. The number c(v) measures 
the extra cost of transporting a mass of m (u) on the system G from the source O 
to the vertex u via the vertex v. 

substep 4: Find the maximum of c (v) over all v € L (u) \ {p(u)}. If max c (v) > 
a [m (u)] , then we find a better parent for the vertex u. In this case, suppose the 
maximum of c(v) is achieved at v* . Then, let 

G* = R '.(G;—m (u) ,v*J + m (u) [v* , u] . 

That is, we change the parent of u from p(u) to v* and then add a mass m (u) at the 
vertex u to the modified path. For convenience, we still denote the final modified 
transport path G* by G. 

step 5: Repeat steps 2-4 until one can not reduce the cost any further. 

4. Examples 

Example 4.1. Let {yi\ be 50 random points in the square [0, 1] x [0, 1]. Then, {y^} 

determines an atomic probability measure 

50 -. 

2 — 1 

Let a= So where O = (0,0) is the origin. Then an optimal transport path from a 
to b looks like the following figures with a — 1,0.75,0.5 and 0.25 respectively: 
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ct=1 a=0.75 




0.2 0.4 0.6 0.8 1 



0.2 0.4 0.6 0.8 1 



Example 4.2. Let {y^} be 100 random points in the rectangle [—2.5,2.5] x [0,1]. 
Then, {yi\ determines an atomic probability measure b = y^._ 1 jkrjS yi . Let a = 8q 
where O = (0,0) is the origin, and let a — 0.85. Then an optimal transport path 
from a to b looks like the following figure. 



1.5- 




-2 -1.5 -1 -0.5 0.5 1 1.5 2 
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Example 4.3. Optimal transport paths from the center to the unit circle. Here, 

the unit circle is represented by 400 points uniformly distributed on the circle. The 

parameter a = 0.75 in the first figure and a = 0.95 in the second one. 

400 points, alpha=0.75, cost=2.0390 400 points, alpha=0.95, cost=1.1989 



0.5 



-0.5 




0.5 



-0.5 




Example 4.4. Optimal transport paths from the center to the unit disk. The first 
one is using random generated points in the disk with a = 2/3 while the second one 
use uniformly generated points in the disk with a — 0.75. 





-1 -0.5 0.5 1 



Example 4.5. An optimal transport path from a point on the boundary to the unit 

square, which is represented by 400 randomly generated points, with a = 0.85. 
a equals b. 3b 




0.2 0.4 0.6 0.8 1 
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Example 4.6. An optimal transport path modeling blood vessels in a placenta of 
new baby 

optimal transport path for Placental ID 1713 
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